%Model Solution with Loan Portfolio, i.e., two loans

%Solve model given q and iterate over q
qSpace=linspace(0,qBar,qGrid);

for z=1:qGrid
        

        %Calculate solution after first loan defaults, which is the
        %baseline solution
        q=qSpace(z);
        SB;
        dFStart=0;
        if (aB<0.0001)
        
                dFStart=(FB-phi*(gamma-r))/phi;
        
        end

        
        solPost = ode15s(@ffunODE,[VB max(kappa*q,VB+0.001)+0.001*(kappa*q==VB)],[FB dFStart-tolerance],options, r, gamma ,Lambda, phi,q,aBar,delta);
        
        %Calculate solution prior to first loan default, taking solution
        %post-default as given
        xSpace=linspace(min(solPost.x),max(solPost.x),1000);
        vFPostEval=deval(solPost,xSpace);
        SBPortfolio;
        
        
        
        
        sol = ode15s(@ffunPortfolio,[VBPre max(kappa*q,VBPre+0.001)+0.001*(kappa*q==VBPre)],[FBPre dFStart-tolerance],options,r, gamma ,Lambda, phi,q,aBar,delta,VB,xSpace,vFPostEval(2,:),vFPostEval(1,:),FB);
        
        %Save initial value in vector
        [m,ind]=min(abs(sol.x-kappa*q));
        vF0(z)=sol.y(1,ind)-0.5*kappa*q^2;

        if (vF0(z)<1.5)

            break;

        end

end

%Determine optimal q
[m,ind]=max(vF0);
q=qSpace(ind);

%Recalculate solution under optimal q
SB;
dFStart=0;

if (aB<0.0001)
        
         dFStart=(FB-phi*(gamma-r))/phi;
        
end
        
solPost = ode15s(@ffunODE,[VB max(kappa*q,VB+0.001)+0.001*(kappa*q==VB)],[FB dFStart-tolerance],options, r, gamma ,Lambda, phi,q,aBar,delta);
        
        
SBPortfolio;
        
        
        

        
        
%Evalvulate value function and calculate payoff-relevant model quantities
xSpace=linspace(min(solPost.x),max(solPost.x),10000);
vFPostEval=deval(solPost,xSpace);
sol = ode15s(@ffunPortfolio,[VBPre max(kappa*q,VBPre+0.001)+0.001*(kappa*q==VBPre)],[FBPre dFStart-tolerance],options,r, gamma ,Lambda, phi,q,aBar,delta,VB,xSpace,vFPostEval(2,:),vFPostEval(1,:),FB);
        


[m,ind]=min(abs(sol.x-kappa*q));



F=sol.y(1,:);
dF=sol.y(2,:);

%Calculate model quantities
for i=1:length(sol.x)
    
    V=sol.x(i);
    F1=sol.y(1,i);
    dF1=sol.y(2,i);
    [m,ind]=min(abs(2*solPost.y(2,:)-dF1));

        

    VPost=solPost.x(ind);
    FPost=solPost.y(1,ind);

    [m,ind]=min(abs(2*vFPostEval(2,:)-dF1));
    VPost=xSpace(ind);
    FPost=vFPostEval(1,ind);



        
        a=(2*(F1-FPost)-(gamma-r)*phi-dF1*(2*V+phi)+dF1*VPost)/(2*phi);

        if (a<0)
            
            a=0;
          
            
        end
        
        if (a>aBar)
            
            a=aBar;
           
            
        end
        
        lambda=Lambda-a-q;
        W=phi*a;
        dV=(gamma+2*lambda)*V-W-lambda*VPost;
        dF2=-(gamma-r)*dF1/(min(dV,-0.0001));

        dW=(-2*V*dF2-dF2*phi+dF2*VPost)/2;
        dWt=dW*dV;
        c=(gamma+2*lambda)*W+phi*a^2-dWt;

        
        vA(i)=a;
        vVPost(i)=VPost;
        vm(i)=m;
        vdV(i)=dV;
        vW(i)=W;
        vc(i)=c;
        vlambda(i)=lambda;
        
end

a0=vA(end);

